{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Parallel and Distributed Optimization\n",
    "\n",
    "This notebook provides concise BlueFog demos to concepts or algorithms introduced in the tutorial. \n",
    "\n",
    "\n",
    "### 1.1 Finite-sum optimizaiton example: distributed least square\n",
    "\n",
    "Suppose $m$ computing nodes collaborate to solve the following problem:\n",
    "\n",
    "\\begin{align*}\n",
    "\\min_{x\\in \\mathbb{R}^d}\\ \\sum_{i=1}^n h_i(x) \\quad \\mbox{where} \\quad h_i(x) = \\frac{1}{2}\\|A_i x - b_i\\|^2 \\hspace{1cm} \\mbox{(Opt-Problem)}\n",
    "\\end{align*}\n",
    "\n",
    "where $h_i(x): \\mathbb{R}^d \\to \\mathbb{R}$ is a local cost function held by node $i$ and $\\{A_i, b_i\\}$ are local data. Each node $i$ can evaluate its own data and gradient, but it has to communicate to achieve information from the other node. We let $x^\\star$ denote the global solution to the above problem\n",
    "\n",
    "### 1.2 All-reduce\n",
    "\n",
    "BlueFog supports the Ring-Allreduce operation as follows: \n",
    "```\n",
    "x_bar = bf.allreduce(x_local)\n",
    "```\n",
    "\n",
    "In the following code, we activate 8 computing nodes (CPUs) and label them as $0,1,\\cdots,7$. We target to let all nodes collaborate to compute the global average of their labels $(\\sum_{i=0}^7 i)/8 = 3.5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-30T05:20:53.953031Z",
     "start_time": "2021-03-30T05:20:53.949476Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting Allreduce.py\n"
     ]
    }
   ],
   "source": [
    "%%writefile Allreduce.py\n",
    " \n",
    "import bluefog.torch as bf\n",
    "import torch\n",
    "\n",
    "bf.init()\n",
    "x_local = torch.ones(1) * bf.rank()\n",
    "x_bar = bf.allreduce(x_local)\n",
    "print('Node {} achieved the global average {}'.format(bf.rank(), x_bar[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "8"
      ]
     },
     "metadata": {},
     "execution_count": 3
    }
   ],
   "source": [
    "# You can change NUM_PROC into any value smaller than the number of CPUs you have.\n",
    "import os\n",
    "NUM_PROC = 2 if os.getenv(\"TEST_ENV\") else 8\n",
    "NUM_PROC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-30T05:20:55.825361Z",
     "start_time": "2021-03-30T05:20:54.582649Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Node 0 achieved the global average 3.5\r\n",
      "Node 1 achieved the global average 3.5\r\n",
      "Node 2 achieved the global average 3.5\r\n",
      "Node 4 achieved the global average 3.5\r\n",
      "Node 5 achieved the global average 3.5\r\n",
      "Node 6 achieved the global average 3.5\r\n",
      "Node 7 achieved the global average 3.5\r\n",
      "Node 3 achieved the global average 3.5\r\n"
     ]
    }
   ],
   "source": [
    "! bfrun -np $NUM_PROC python Allreduce.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Distributed Gradient Descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The distributed gradient descent is $x^{k+1} = x^k - \\frac{\\alpha}{n}\\sum_{i=1}^n \\nabla h_i(x^k)$. This process can be implemented as\n",
    "\n",
    "```python\n",
    "# core distributed gradient descent snippet \n",
    "grad = A.T.mm(A.mm(x) - b)\n",
    "next_x = x - alpha * bf.allreduce(grad)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-30T05:20:56.918844Z",
     "start_time": "2021-03-30T05:20:56.913919Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting DistributedGD.py\n"
     ]
    }
   ],
   "source": [
    "%%writefile DistributedGD.py\n",
    " \n",
    "import bluefog.torch as bf\n",
    "import torch\n",
    "\n",
    "bf.init()\n",
    "# Make sure different agent has different random seed.\n",
    "torch.manual_seed(12345 * bf.rank())\n",
    "\n",
    "def generate_data(m, d):\n",
    "    A = torch.randn(m, d).to(torch.double)\n",
    "    ns = 0.1*torch.randn(m, 1).to(torch.double)\n",
    "    x_o = torch.rand(d,1).to(torch.double)\n",
    "    b = A.mm(x_o) + ns\n",
    "    \n",
    "    return A, b\n",
    "\n",
    "def check_opt_cond(x, A, b):\n",
    "    \n",
    "    grad_local = A.t().mm(A.mm(x) - b)\n",
    "    grad = bf.allreduce(grad_local, name='gradient')  # global gradient\n",
    "    \n",
    "    # the norm of global gradient is expected to be 0 (optimality condition)\n",
    "    global_grad_norm = torch.norm(grad, p=2)\n",
    "    print(\"[Distributed Grad Descent] Rank {}: global gradient norm: {}\".format(bf.rank(), global_grad_norm))\n",
    "        \n",
    "    return\n",
    "    \n",
    "\n",
    "def distributed_grad_descent(A, b, maxite=5000, alpha=1e-1):\n",
    "\n",
    "    m, d = A.shape\n",
    "    \n",
    "    x_opt = torch.zeros(d, 1, dtype=torch.double)\n",
    "\n",
    "    for _ in range(maxite):\n",
    "        # calculate local gradient \n",
    "        grad_local = A.t().mm(A.mm(x_opt) - b)\n",
    "        \n",
    "        # global gradient\n",
    "        grad = bf.allreduce(grad_local, name='gradient')\n",
    "\n",
    "        # distributed gradient descent\n",
    "        x_opt = x_opt - alpha*grad\n",
    "\n",
    "    check_opt_cond(x_opt, A, b)\n",
    "    \n",
    "    return x_opt\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    m, d = 20, 5 # dimension of A\n",
    "    A, b = generate_data(m, d)\n",
    "    x_opt = distributed_grad_descent(A, b, maxite=200, alpha=1e-2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $x^k\\to x^\\star$, it holds that $\\sum_{i=1}^m \\nabla h_i(x^k) \\to 0$. We can use $\\|\\sum_{i=1}^m \\nabla h_i(x)\\|$ as a metric to gauge whether distributed gradient descent converge or not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-30T05:20:59.456793Z",
     "start_time": "2021-03-30T05:20:57.890226Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Distributed Grad Descent] Rank 2: global gradient norm: 7.786811775375865e-15\r\n",
      "[Distributed Grad Descent] Rank 3: global gradient norm: 7.786811775375865e-15\r\n",
      "[Distributed Grad Descent] Rank 4: global gradient norm: 7.786811775375865e-15\r\n",
      "[Distributed Grad Descent] Rank 5: global gradient norm: 7.786811775375865e-15\r\n",
      "[Distributed Grad Descent] Rank 6: global gradient norm: 7.786811775375865e-15\r\n",
      "[Distributed Grad Descent] Rank 7: global gradient norm: 7.786811775375865e-15\r\n",
      "[Distributed Grad Descent] Rank 0: global gradient norm: 7.786811775375865e-15\r\n",
      "[Distributed Grad Descent] Rank 1: global gradient norm: 7.786811775375865e-15\r\n"
     ]
    }
   ],
   "source": [
    "! bfrun -np $NUM_PROC python DistributedGD.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4 Primal Decomposition\n",
    "\n",
    "We consider the following linearly-constrained resource sharing problem:\n",
    "\n",
    "\\begin{align*}\n",
    "\\min_{\\{x_i\\}, y}&\\quad \\|y\\|_1 + \\frac{1}{2}\\sum_{i=1}^n \\|A_i x_i - b_i\\|^2 \\\\\n",
    "\\mathrm{s.t.} &\\quad\\sum_{i=1}^n B_i x_i = y\n",
    "\\end{align*}\n",
    "\n",
    "where $x_i \\in \\mathbb{R}^d$, $A_i \\in \\mathbb{R}^{m\\times d}$, and $B_i \\in \\mathbb{R}^{d\\times d}$. For simplicity, we assume each $B_i$ is a invertible matrix. Each node $i$ has local data $\\{A_i, B_i, b_i\\}$. To solve the above problem in a distributed manner, we introduce $y_i = B_i x_i$ so that $y=\\sum_{i=1}^n y_i$. Since $B_i$ is invertible, we have $x_i = B_i^{-1} y_i$. Substituting these facts into the above problem, we achieve\n",
    "\n",
    "\\begin{align*}\n",
    "\\min_{\\{y_i\\}}&\\quad \\|y_1 + \\cdots + y_n\\|_1 + \\frac{1}{2}\\sum_{i=1}^n \\|A_i B_i^{-1} y_i - b_i\\|^2 \n",
    "\\end{align*}\n",
    "\n",
    "which can be solved in the following distributed manner. For notation simplicity, we let $C_i = A_i B_i^{-1}$ and $g(y_1,\\cdots, y_n) = \\|y_1 + \\cdots + y_n\\|_1$. Following proximal gradient descent, we have \n",
    "\n",
    "\\begin{align*}\n",
    "z_i^{k+1} &= y_i^k - \\alpha C_i^T(C_i y_i^k - b_i), \\quad \\forall\\ i=1,\\cdots, n\\\\\n",
    "y_i^{k+1} &= [\\mathrm{Prox}_{\\alpha g}(z_1^{k+1},\\cdots, z_n^{k+1})]_{i} =z_i^{k+1} - \\alpha v^{k+1}, \\quad \\forall\\ i=1,\\cdots, n \\\\\n",
    "x_i^{k+1} &= B_i^{-1} y_i^{k+1}, \\quad \\forall\\ i=1,\\cdots, n\n",
    "\\end{align*}\n",
    "\n",
    "where \n",
    "$$v^{k+1} = \\frac{1}{n\\alpha^2}\\big(\\alpha(z^{k+1}_1 + \\cdots + z^{k+1}_n) - \\mathrm{Prox}_{n \\alpha^2 \\|\\cdot\\|_1}(\\alpha z^{k+1}_1 + \\cdots + \\alpha z^{k+1}_n)\\big) \\quad \\mbox{(Need Allreduce Communication)}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-30T05:20:59.628432Z",
     "start_time": "2021-03-30T05:20:59.622147Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting PrimalDecompose.py\n"
     ]
    }
   ],
   "source": [
    "%%writefile PrimalDecompose.py\n",
    " \n",
    "import bluefog.torch as bf\n",
    "import torch\n",
    "\n",
    "bf.init()\n",
    "# Make sure different agent has different random seed.\n",
    "torch.manual_seed(12345 * bf.rank())\n",
    "\n",
    "def generate_data(m, d):\n",
    "    A = torch.randn(m, d).to(torch.double)\n",
    "    B_inv = torch.randn(d, d).to(torch.double)\n",
    "    ns = 0.1*torch.randn(m, 1).to(torch.double)\n",
    "    x_o = torch.rand(d,1).to(torch.double)\n",
    "    b = A.mm(x_o) + ns\n",
    "    \n",
    "    return A, B_inv, b\n",
    "\n",
    "def soft_threshold(x, kappa):\n",
    "    \n",
    "    zeros = torch.zeros(d,1).to(torch.double)\n",
    "    x = torch.max(x - kappa, zeros)\n",
    "    x = torch.min(x + kappa, zeros)\n",
    "    \n",
    "    return x\n",
    "\n",
    "def check_opt_cond(y, A, B_inv, b, alpha):\n",
    "    \n",
    "    m, d = A.shape\n",
    "    n = bf.size()\n",
    "    C = A.mm(B_inv)\n",
    "    \n",
    "    # update z\n",
    "    grad_local = C.t().mm(C.mm(y) - b)\n",
    "    z = y - alpha*grad_local\n",
    "    z_bar = bf.allreduce(z)\n",
    "    z_sum = z_bar * n\n",
    "\n",
    "    # update v\n",
    "    v = (1/(n*alpha*alpha)) * (alpha*z_sum - soft_threshold(alpha*z_sum, n*alpha*alpha))\n",
    "\n",
    "    # update y\n",
    "    y_next = z - alpha*v\n",
    "    \n",
    "    # the norm of global gradient is expected to be 0 (optimality condition)\n",
    "    global_grad_norm = torch.norm((y - y_next), p=2)\n",
    "    print(\"[Primal Decomposition] Rank {}: optimality metric norm: {}\".format(bf.rank(), global_grad_norm))\n",
    "        \n",
    "    return\n",
    "\n",
    "def primal_decomposition(A, B_inv, b, maxite=5000, alpha=1e-1):\n",
    "\n",
    "    m, d = A.shape\n",
    "    n = bf.size()\n",
    "    C = A.mm(B_inv)\n",
    "    \n",
    "    y = torch.zeros(d, 1, dtype=torch.double)\n",
    "\n",
    "    for _ in range(maxite):\n",
    "        # update z\n",
    "        grad_local = C.t().mm(C.mm(y) - b)\n",
    "        z = y - alpha*grad_local\n",
    "        z_bar = bf.allreduce(z)\n",
    "        z_sum = z_bar * n\n",
    "        \n",
    "        # update v\n",
    "        v = (1/(n*alpha*alpha)) * (alpha*z_sum - soft_threshold(alpha*z_sum, n*alpha*alpha))\n",
    "        \n",
    "        # update y\n",
    "        y = z - alpha*v\n",
    "\n",
    "        # update x\n",
    "        x = B_inv.mm(y)\n",
    "\n",
    "    check_opt_cond(y, A, B_inv, b, alpha)\n",
    "    \n",
    "    return y\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    m, d = 20, 5 # dimension of A\n",
    "    maxite = 5000 if os.getenv(\"TEST_ENV\") else 100\n",
    "    A, B_inv, b = generate_data(m, d)\n",
    "    y = primal_decomposition(A, B_inv, b, maxite=maxite, alpha=3e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define $G^k_i = \\frac{1}{\\alpha}\\big(y_i^k - [\\mathrm{Prox}_{\\alpha g}(z_1^{k+1},\\cdots, z_n^{k+1})]_i\\big)$ and use $\\|G^k_i\\|$ as the metric to evaluate the optimality of $y_i^k$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-30T05:21:15.021227Z",
     "start_time": "2021-03-30T05:21:00.552256Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Primal Decomposition] Rank 2: optimality metric norm: 1.526888906516414e-05\r\n",
      "[Primal Decomposition] Rank 3: optimality metric norm: 0.00010111695583238024\r\n",
      "[Primal Decomposition] Rank 5: optimality metric norm: 2.66226578246676e-06\r\n",
      "[Primal Decomposition] Rank 6: optimality metric norm: 3.1143558672548934e-05\r\n",
      "[Primal Decomposition] Rank 4: optimality metric norm: 1.2605012098981551e-05\r\n",
      "[Primal Decomposition] Rank 1: optimality metric norm: 0.00013532193550425833\r\n",
      "[Primal Decomposition] Rank 0: optimality metric norm: 6.13283876701975e-05\r\n",
      "[Primal Decomposition] Rank 7: optimality metric norm: 4.335475665174773e-06\r\n"
     ]
    }
   ],
   "source": [
    "! bfrun -np $NUM_PROC python PrimalDecompose.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.5 Distributed ADMM\n",
    "\n",
    "Suppose we want to optimize the following consensus problem:\n",
    "$$\n",
    "    \\min_{x} \\sum_{i=1}^n h_i(x)\n",
    "$$\n",
    "Each machine $i$ can access the function $h_i$ only. Note the problem has a shared $x$ cross all machines. To relax this condition, we can reformat the problem into the following equivalent form:\n",
    "\\begin{align}\n",
    "    \\min_{\\{x_i\\}, y} & \\;\\;\\; \\sum_{i=1}^n h_i(x_i),\\\\\n",
    "    {\\rm subject\\ to} & \\;\\;\\;  x_i = y,\\;\\; \\forall i. \n",
    "\\end{align}\n",
    "\n",
    "Applying ADMM, we obtain:\n",
    "\n",
    "\\begin{align}\n",
    "    x^{k+1}_i =&\\;\\; {\\rm argmin}_{x_{i}} \\left\\{ h_i(x_i) + \\langle u_i^k, x_i - y^k\\rangle + \\frac{\\alpha}{2}  \\left\\|x_i - y^k \\right\\|^2\\right\\} \\\\\n",
    "    y^{k+1} =&\\;\\; \\frac{1}{n} \\left(x_i^{k+1} + \\frac{1}{\\alpha} u^k_i\\right)\\\\\n",
    "    u^{k+1}_i = &\\;\\; u^k_i + \\alpha\\left(x_i^{k+1} - y^{k+1} \\right)\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can further simplify it by noticing that $u^k_1,\\ldots,u^k_n$ has mean 0. Here we provide the code snippet of a simple quadratic cost function again. \n",
    "\n",
    "``` python\n",
    "# For each agent, it owns data A and b differently. \n",
    "def CentralizedADMMStepL2(A, b, x, y, u, alpha):\n",
    "    next_x = ProximalStepL2(A, b, x, y, u, alpha)\n",
    "    # We use allreduce to mimic the centralized behavior\n",
    "    # It should be based on PS architecture and using gather and broadcast.\n",
    "    next_y = bf.allreduce(next_x)   # Without u is okay since allreudce(u) == 0\n",
    "    next_u = u + alpha * (next_x - next_y)\n",
    "    return next_x, next_y, next_u\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-03-30T05:21:17.723147Z",
     "start_time": "2021-03-30T05:21:15.026908Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Centralized ADMM] Rank 5: ADMM residue gradient norm: 5.6263947243994325e-15\r\n",
      "[Centralized ADMM] Rank 0: ADMM residue gradient norm: 5.6263947243994325e-15\r\n",
      "[Centralized ADMM] Rank 1: ADMM residue gradient norm: 5.6263947243994325e-15\r\n",
      "[Centralized ADMM] Rank 2: ADMM residue gradient norm: 5.6263947243994325e-15\r\n",
      "[Centralized ADMM] Rank 4: ADMM residue gradient norm: 5.6263947243994325e-15\r\n",
      "[Centralized ADMM] Rank 6: ADMM residue gradient norm: 5.6263947243994325e-15\r\n",
      "[Centralized ADMM] Rank 3: ADMM residue gradient norm: 5.6263947243994325e-15\r\n",
      "[Centralized ADMM] Rank 7: ADMM residue gradient norm: 5.6263947243994325e-15\r\n",
      "Last three entries of x_ar:\r\n",
      " tensor([[0.4883],\r\n",
      "        [0.6776],\r\n",
      "        [0.6447]], dtype=torch.float64)\r\n",
      "Last three entries of x_admm:\r\n",
      " tensor([[0.4883],\r\n",
      "        [0.6776],\r\n",
      "        [0.6447]], dtype=torch.float64)\r\n"
     ]
    }
   ],
   "source": [
    "! bfrun -np $NUM_PROC python CentralizedADMM.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEWCAYAAACe8xtsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hUVfrA8e+bEAgQeglKMaETOgm9BWkBAUFQiiKySFtREXWV1Z/iumtZG6BYQBArAZG+CNJCkR6QJl0QAghKk9AJ7++PGdiQTZkMSSYzeT/Pcx9yz9zyHi7knXvOveeIqmKMMcakh5+nAzDGGON9LHkYY4xJN0sexhhj0s2ShzHGmHSz5GGMMSbdLHkYY4xJN0sexmQQEYkUkbhE6ztEJDKDzzFZRP6Zkcc0xh2WPIxXEZE+IrJRROJF5JiIfC8izTLguKNE5KuMiPEGVa2uqjEZeUxXOJOYishzScpDnOXxzuW4iMwTkbZJtjsoIldEpHiS8s3O/UOc65Od6/cm2e49Z/kjmVJBky1Y8jBeQ0RGAKOB14BgoBzwIXBvavtl0LlFRLzl/0s/4BTwcAqfF1bVIKA2sAiYmcwv+gNA7xsrIlITyJfMsfYkPo+I5AIeAPa7G7zxDt7yn8HkcCJSCPgH8JiqzlDV86p6VVXnquqzzm38ROR5EdkvIidFZJqIFHV+duNbdz8ROSQif4jIC87PooC/Az2d38i3OMtjRORfIvIjcAEoLyL9RWSniJwTkV9EZHAqMR8UkTbOn88k+sZ/Psk3+E4i8pNzm9UiUivRMeqKyCbn+aYCgWn8PeUHegCPAZVEJCKlbVX1N1UdA4wC3kySHL/k1uTTD/gimcPMBZqJSBHnehSwFfgttTiN97PkYbxFYxy/OGemss3jQFegJXAncBoYl2SbZkAVoDXwkohUU9UFOO5mpqpqkKrWTrR9X2AQUAD4FTgBdAIKAv2B90SkXlrBq2ph57GDgDHASuCIiNQFJgGDgWLAJ8AcEckjIrmBWTh+kRcFvgW6p3Gq+4B457YLcfzST8sMoCSOv5cb1gIFRaSaiPgDvYDkmvUuAbOdn4Mj4SSXZIyPseRhvEUx4A9VvZbKNkOAF1Q1TlUv4/hG3cPZlHLDK6p6UVW3AFtwNN2kZrKq7lDVa847nf+o6n51WA78ADR3tRIi0hPoA3RX1as4EtMnqrpOVRNU9XPgMtDIuQQAo53nng5sSOMU/XAkwQTgG6CXiASksc9R559Fk5TfuPtoC+wEjqSw/xfAwyJSGEfinpXG+YwPsORhvMVJoHiSRJDUXTja78+IyBkcv/AScPSP3JC4OeUCEJTGeQ8nXhGRDiKyVkROOc/RESie/K63ct5lfAB0U9XfE8X89I2Ynccsi+PO6U7giN46eumvqRy/LNAK+NpZNBvH3do9aYRW2vnnqSTlX+JIdI+Qyt2Eqq4CSgAvAPNU9WIa5zM+wJKH8RZrcHwj75rKNoeBDs4mohtLoKqm9I05sZSGl75ZLiJ5gO+At4FgVS0MzAckrYOLSEkc38gfU9XNSWL+V5KY86nqFOAYUFpEEh+/XCqn6Yvj//RcEfkN+AVH8kir6aobjua43YkLVfVXHB3nHXE0baXmK+BprMkqx7DkYbyCqp4FXgLGiUhXEcknIgHOO4F/Ozf7GPiXiNwFICIlkj5GmorjQEgaT1TlBvIAvwPXRKQD0C6tAzvvlqYDX6nqtCQfTwCGiEhD5xNd+UXkHhEpgCNhXgOecNb1PqBBKqfqB7wC1Em0dAc6ikixZOIKFpFhwMvASFW9nswxBwB3q+r5NKo5Fkfz1oo0tjM+wpKH8Rqq+g4wAngRxy/ww8Aw/tvGPgaYA/wgIudwdPo2dPHw3zr/PCkim1I4/zngCWAajs74Ps7zpaUMjn6R4YmeuIoXkXKquhEYiKM56zSwD0czEap6BUcH+CM4mpR6ksIdgIg0wtEENs75FNWNZY7zmL0TbX5GRM4D23DcVdyvqpNSqPN+Z4ypUtVTqrokSROb8WFi19oYY0x62Z2HMcaYdLPkYYwxJt0seRhjjEk3Sx7GGGPSLbUXrnxG8eLFNSQkxO39z58/T/78+TMuoGwup9UXrM45hdU5fWJjY/9Q1RLJfZYjkkdISAgbN6b5tGGKYmJiiIyMzLiAsrmcVl+wOucUVuf0EZEURzSwZitjjDHpZsnDGGNMulnyMMYYk245os/DGJO2q1evEhcXx6VLlzwdSqYpVKgQO3fu9HQYWcqVOgcGBlKmTBkCAtIavf+/LHkYYwCIi4ujQIEChISEcOtAvr7j3LlzFChQwNNhZKm06qyqnDx5kri4OEJDQ10+rjVbGWMAuHTpEsWKFfPZxGGSJyIUK1Ys3XecljyMMTdZ4siZ3LnuljxSceXadV6Zu4Mzl5Ob5sAYY3IuSx6p2HnsT6asP8Q/1lxi+5Gzng7HGJ/322+/0atXLypUqEB4eDgdO3Zkz549bh1r8uTJHD16NO0Nkxg1ahRvv/02AC+99BKLFy926/yJBQWlPNvxrFmzEBF27dp1s+zgwYPkzZuXunXrUq1aNRo0aMDkyZNvfj558mRE5JbYbhxn+vTpAERGRlKuXDkST7vRtWvXVGNJD0seqahdtjDThzQBoMfHq/nP1mMejsgY36WqdOvWjcjISPbv309sbCyvv/46x48fd+t4qSWPhIQEl47xj3/8gzZt2rh1fldNmTKFZs2aMWXKlFvKK1SowObNm9m5cyfR0dGMHj2azz777ObnNWvWJDo6+pbj1K5d+5ZjFC5cmLVr1wJw5swZjh3LuN9hljzSUKN0IV5unJfqdxbisW82MXXDIU+HZIxPWrZsGQEBAQwZMuRmWe3atWnevDkAb731FvXr16dWrVq8/PLLgOMberVq1Rg4cCDVq1enXbt2XLx4kenTp7Nx40YefPBB6tSpw8WLFwkJCeGll16iXr16fPvtt0yYMIH69etTu3ZtunfvzoULF/4npkceeeTmserUqUOdOnWoWbPmzT6C/fv3ExUVRXh4OM2bN79593DgwAEaN25MzZo1efHFF1Osc3x8PKtWrWLixIm3JIKkypcvz7vvvsvYsWNvljVv3pz169dz9epV4uPj2bdvH3Xq1Lllv169et28E5kxYwb33XdfqtcgPexRXRcUyiN8M7AhA7+I5cVZ2ylfIoj6IUU9HZYxmeaVuTv4+eifGXrMsDsL8nLn6il+vn37dsLDw5P97IcffmDv3r2sX78eVaVLly6sWLGCcuXKsXfvXqZMmcKECRN44IEH+O6773jooYf44IMPePvtt4mIiLh5nKJFi7Jpk2OW4ZMnTzJw4EAAXnzxRSZOnMjjjz+e7PkjIiL46aefAHj22WeJiooCYNCgQXz88cdUqlSJdevW8de//pWlS5fy5JNPMnToUB5++GHGjRuXYp1nz55NVFQUlStXplixYsTGxqb4d1CvXr1bmrZEhDZt2rBw4ULOnj1Lly5dOHDgwC37tG7dmgEDBpCQkEB0dDTjx4/n1VdfTTGe9LA7DxflyeXP+73qUqZIPoZ+FcuRMxc9HZIxOcYPP/zADz/8QN26dW/+Et27dy8AoaGhN79xh4eHc/DgwRSPk/ib9/bt22nevDk1a9bk66+/ZseOHWnGMXXqVDZt2sQbb7xBfHw8q1ev5v7776dOnToMHjz4ZrPQjz/+SO/ejmnj+/btm+LxpkyZQq9evQDHXULSpqvEkpsyvFevXkRHRxMdHX3zfIn5+/vTqFEjoqOjb959ZRSvvPMQkfzAh8AVIEZVv86K8xbKF8CEh8PpOm41g77YyPQhTcib2z8rTm1MlkrtDiGzVK9e/WYTS1KqysiRIxk8ePAt5QcPHiRPnjw31/39/bl4MeUvdomHJn/kkUeYNWsWtWvXZvLkycTExKQa3/bt2xk1ahQrVqzA39+f69evU7hw4Zt3JEml9fjrqVOnWLp0Kdu2bUNESEhIQER46623kt1+8+bNVKtW7ZayBg0asG3bNvLly0flypWT3a9Hjx48+OCDjBo1KtV40ivb3HmIyCQROSEi25OUR4nIbhHZJyLPO4vvA6ar6kCgS1bGWbFkAcb2rsPPx/5kxLSfuH79f78NGGPS7+677+by5cuMHz/+ZtnWrVtZuXIl7du3Z9KkScTHxwNw5MgRTpw4kerxChQowLlz51L8/Ny5c9xxxx1cvXqVr79O/fvnmTNn6N27N1988QUlSjimtyhYsCChoaF8++23gCPBbdmyBYCmTZve7MNI6djTp0+nb9++/Prrrxw8eJDDhw8TGhrKypUr/2fbgwcP8swzzyTbrPbGG2/w2muvpRh7kyZNGDlyZLJ3Jrcj2yQPYDIQlbhARPyBcUAHIAzoLSJhQBngsHMz1x6byEB3Vw3mhY7V+H77b7yxYFfaOxhj0iQizJw5k8WLF1OhQgWqV6/OyJEjKVWqFO3ataNPnz43O6F79OiRamIAx53FkCFDbnaYJ/Xqq6/SsGFDmjZtStWqVVM91uzZs/n1118ZOHDgzY5zcCSGiRMnUrt2bapXr87s2bMBGDNmDOPGjaNmzZocOXIk2WNOmTKFbt263VLWvXv3m01X+/fvv/mo7gMPPMATTzxB//79/+c4HTp0oFWrVinGLiI888wzFC9ePNU6ppck147mKSISAsxT1RrO9cbAKFVt71wf6dw0DjitqvNEJFpVeyVzrEHAIIDg4ODw1J5kSEt8fPz/PButqny18wpLDl2jb1huWpdzfUCx7C65+vo6q7NjAL2KFSt6MKLMl5CQgL9/zmpqdrXO+/bt4+zZW99na9WqVayqRiS3fXbv8yjNf+8wwJE0GgJjgQ9E5B5gbnI7qup4YDxARESE3s7sYSnNxNWipTLoi418vfMELSNq0SYs2O1zZCc221rOkLTOO3fu9PlBA21gxJQFBgZSt25dl4+bnZqtXKaq51W1v6oOzarO8uT4+wlje9elRmnHOyDrfjnpqVCMMSZLZffkcQQom2i9jLMs28ifJxeT+zegTJG8PPr5RhvGxHi17NSMbbKOO9c9uyePDUAlEQkVkdxAL2COh2P6H0Xz5+arRxtSMG8A/SatZ9+JeE+HZEy6BQYGcvLkSUsgOcyN+TwCAwPTtV+26fMQkSlAJFBcROKAl1V1oogMAxYC/sAkVU37TR4PuKNQXr4c0IAHPllDr/FrmTKwIZWCc1bbqvFuZcqUIS4ujt9//93ToWSaS5cupfuXpLdzpc43ZhJMj2yTPFQ12YeQVXU+MD+Lw3FL+RJBTBnYiD6frqPX+LV89WhDqt1R0NNhGeOSgICAdM0k541iYmLS1SnsCzKrztm92crrVAouwLTBjcmdy4/eE9ayNe6Mp0MyxpgMZ8kjE4QWz8+0wY0JypOLXuPXsmxX6m/CGmOMt7HkkUnKFs3HjKFNKF8iP49+sZHo9TaUuzHGd1jyyEQlCwYydVBjmlUszvMztvHvBbtsLCxjjE+w5JHJ8ufJxaf9IujdoBwfxuxn8FexxF++5umwjDHmtljyyAIB/n681q0GozqHsXTXCbp/uJrDp/531jJjjPEWljyyiIjwSNNQJvevz7GzF+n8wSpW7f3D02EZY4xbLHlkseaVSjBnWDNKFsjDw5PWMX7Ffnuj1xjjdSx5eEBI8fzM/GtTomqU4rX5uxg2ZTPnrR/EGONFLHl4SP48uRjXpx7PRVXl+23H6Pbhjxz447ynwzLGGJdY8vAgEWFoZAW++EtDfj93mS7vr2Lxz8c9HZYxxqTJkkc20KxSceY+3oyQ4o4XCt9euJsEex/EGJONWfLIJsoUyce3QxrTM6IsHyzbxyOfref0+SueDssYY5JlySMbCQzw580etXjjvpqs++UUnd5fxZbDNrCiMSb7seSRDfVqUI7pQxsDcP/Ha/hm3SF7nNcYk61Y8simapUpzLzHm9GoQjH+PnMbz3y7lYtXEjwdljHGAF6aPESkq4hMEJGpItLO0/FkliL5c/PZI/UZ3qYSMzbH2eO8xphsI8uTh4hMEpETIrI9SXmUiOwWkX0i8nxqx1DVWao6EBgC9MzMeD3N308Y3qYyk/s34Lc/L9Hl/VUs2H7M02EZY3I4T9x5TAaiEheIiD8wDugAhAG9RSRMRGqKyLwkS8lEu77o3M/ntaxcgv880ZzyJYMY8tUm/jnvZ64mXPd0WMaYHEo80RErIiHAPFWt4VxvDIxS1fbO9ZEAqvp6CvsL8AawSFUXp7DNIGAQQHBwcHh0dLTb8cbHxxMUFOT2/hnp2nUletcVFh+6RsXCfvy1Th6KBmbsd4DsVN+sYnXOGazO6dOqVatYVY1I7rNctxVVxikNHE60Hgc0TGX7x4E2QCERqaiqHyfdQFXHA+MBIiIiNDIy0u3gYmJiuJ39M1qbu2HulqM8/91W/rkhgfd61qRl5RIZdvzsVt+sYHXOGazOGccrO8xVdayqhqvqkOQSR07QufadzHm8GSWC8vDIZ+t59wd7K90Yk3WyS/I4ApRNtF7GWWZSUaFEELMea0qPemUYu3QffSeu48S5S54OyxiTA2SX5LEBqCQioSKSG+gFzPFwTF4hb25/3rq/Nv/uUYtNh05zz9hVrNl/0tNhGWN8nCce1Z0CrAGqiEiciAxQ1WvAMGAhsBOYpqo7sjo2b/ZARFlmP9aMAoG5ePDTtby/ZC/XrRnLGJNJsrzDXFV7p1A+H5ifxeH4lCqlCjB3WDP+PnMb7yzaw/qDp3ivZx2KB+XxdGjGGB+TXZqtTAbJnycXo3vW4fX7arL+wCk6jllpzVjGmAxnycMHiQi9G5Rj1mNNCcrjaMYas3ivPY1ljMkwljx8WLU7CjLn8WZ0qX0n7y3eY09jGWMyjCUPHxeUJxfv9axz82msjmNWsmLP754Oyxjj5Sx55AAiwgMRZZkzrBlF8+fm4UnreXPBLhsbyxjjNkseOUjl4ALMfqwZvRuU5aOY/fT8ZA2HT13wdFjGGC9kySOHyZvbn9fvq8X7veuy93g8Hceu5D9bbYh3Y0z6WPLIoTrXvpP/PNGcCiWCeOybTYycsZULV655OixjjJew5JGDlSuWj2+HNGZIywpEbzhM5/dXsePoWU+HZYzxApY8crgAfz+e71CVrwY0JP7yNbqNW83Cg1dtaBNjTKoseRgAmlYszoInWxBZpQRTdl2h32frOf6nvRNijEmeJQ9zU5H8ufmkbzj9wnKz4eApokavYOGO3zwdljEmG7LkYW4hIrQqF8C8x5tTukheBn8Zy3PTt3L+snWmG2P+y5KHSVbFkkHMGNqUoZEVmBZ7mI5jVxL762lPh2WMySbSTB4i4i8iy7IiGJO95M7lx3NRVZk6qDHXEpT7P17N2wt3c+WavZluTE6XZvJQ1QTguogUyoJ4TDbUILQoC4Y3p1vdMnywbB/3ffQje4+f83RYxhgPcrXZKh7YJiITRWTsjSUzA0uLiOQXkY0i0smTceQUBQIDeOeB2nz8UDhHz1zinvdX8enKX+yRXmNyKFeTxwzg/4AVQGyiJd1EZJKInBCR7UnKo0Rkt4jsE5HnXTjUc8A0d2Iw7ouqUYqFw1vQolJx/vmfnfSasNbGxzImB3JpGlpV/VxEcgOVnUW7VfWqm+ecDHwAfHGjQET8gXFAWyAO2CAicwB/4PUk+/8FqA38DAS6GYO5DSUK5GHCwxF8GxvHP+b+TNToFbzYKYxe9csiIp4OzxiTBUQ17WYHEYkEPgcOAgKUBfqp6gq3TioSAsxT1RrO9cbAKFVt71wfCaCqSRPHjf3/BeQHwoCLQDdVvZ5km0HAIIDg4ODw6Ohod0IFID4+nqCgILf39zbpqe8fF68zcdtldp66Ts3i/vylRm6KBHrfQ3w57RqD1TmnuJ06t2rVKlZVI5L9UFXTXHA0UVVJtF4ZiHVl3xSOFwJsT7TeA/g00Xpf4AMXjvMI0Cmt7cLDw/V2LFu27Lb29zbprW9CwnX9fPUBrfri91rj5QU6feNhvX79euYEl0ly2jVWtTrnFLdTZ2CjpvB71dWviAGqujtRwtkDBKQvh2U8VZ2sqvM8HUdO5+cnPNw4hO+fbE6V4AI8/e0WBn4Ra1PeGuPDXE0esSLyqYhEOpcJwMYMjOMIjqawG8o4y4wXCSmen6mDG/PiPdVYufd32r23gtk/Hblxl2iM8SGuJo8hODqon3AuPwNDMzCODUAlEQl1dsz3AuZk4PFNFvH3Ex5tXp75TzYntHh+noz+icFf2l2IMb7GpTfMgS2q+q6q3udc3lPVy+6cUESmAGuAKiISJyIDVPUaMAxYCOwEpqnqDneOb7KHCiWCmD6kCX/vWJWYPXYXYoyvSfNRXVVNcL5/UU5VD93uCVW1dwrl84H5t3t8k334+wmDWlTg7qrBPDt9C09G/8TcLcd4rVsNSha0p6yN8WauNlsVAXaIyBIRmXNjyczAjO+oWNJxF/JCR0dfSJt3l/NdbJzdhRjjxVx6SRDH2+XGuM3fTxjYojytq5Xkb9O38vS3W5i39Siv3VeTOwrl9XR4xph0crXP4xNVXZ50yYL4jI8pXyKIaYMb83LnMNb+cop2765gyvpDdhdijJdxdVTd3SJSLgviMTmAn5/Qv2koC4e3oEbpQoycsY2HJq6zMbKM8SLW52E8plyxfHwzsCGvdavJlsNnaffeCiatOkCCjdRrTLZnfR7Go0SEPg3LEVmlBC/M3MY/5v3Mf7Yd483uNalYsoCnwzPGpCDVOw8RqQrg7N9Ym6S/w633PIxJzp2F8zLpkfq817M2+3+Pp+OYVXywdC9XE2zWQmOyo7Sarb5J9POaJJ99mMGxmBxOROhWtwyLnmpJ27Bg3v5hD10++JFtcWc9HZoxJom0koek8HNy68ZkiBIF8jDuwXp80jeck/GX6frhj7z+/U4uXU3wdGjGGKe0koem8HNy68ZkqPbVS7FoREt61CvDJ8t/IWr0CtbsP+npsIwxpN1hXsY5V7kk+hnneulMjcwYoFDeAN7sUYsude5k5Ixt9J6wlt4NyvF8h6oUyuvxWQGMybHSSh7PJvo56RDsGTkkuzGpalqxOAuHt+DdRbuZuOoAS3Ye59WuNWhfvZSnQzMmR0o1eajq50nLRKSUqv6WeSEZk7y8uf154Z4wOtW6k+e+28rgL2PpWLMUo7pUp2QBG2jRmKzkzmTTNvKt8ajaZQsz9/FmPNu+Cot3nqDNO8uZusGGODEmK7mTPOwpK+NxAf5+PNaqIt8/2ZyqdxTkue+20WfCOg78cd7ToRmTI7iTPCZkeBTGuKlCiSCiBzbitW412X70LFGjV/BhzD57udCYTOZy8hCRZiLSX1U/FJESIhKamYGlEYufiPxLRN4XkX6eisNkD35+jiFOloxoyd1VS/LvBbvp8sGPbI074+nQjPFZLiUPEXkZeA4Y6SwKAL5y54QiMklETojI9iTlUc4ZC/eJyPNpHOZeoAxwFYhzJw7je0oWDOSjh8L5+KFwTp2/TNdxP/LqvJ+5cOWap0Mzxue4eufRDegCnAdQ1aOAu6PWTQaiEhc45wwZB3QAwoDeIhImIjVFZF6SpSRQBVitqiOAoW7GYXxUVA3Hy4W9G5Rj4qoDtH13BTG7T3g6LGN8irjyhIqIrFfVBiKySVXriUh+YI2q1nLrpCIhwDxVreFcbwyMUtX2zvWRAKr6egr7PwRcUdVpIjJVVXsms80gYBBAcHBweHR0tDuhAhAfH09QUJDb+3sbX6rvntMJfLb9MsfOK43u8KdP1TwUzPO/z3z4Up1dZXXOGW6nzq1atYpV1YjkPnN1SPZpIvIJUFhEBgJ/IWM7zksDhxOtxwENU9l+BvC+iDQHViS3gaqOB8YDREREaGRkpNvBxcTEcDv7extfqm8k0K9zAh8u28+HMfvYdfYqL3SsRo/wMoj8N4n4Up1dZXXOGTKrzmkmD3H8D5sKVAX+xNFk9JKqLsrwaFykqheAAZ46v/EueXL581TbynSqdQcjZ2zj2elbmbn5CK91q0lI8fyeDs8Yr+TKNLQKzFfVRar6rKo+kwmJ4whQNtF6GWeZMRmmUnABpg1uzD+71mBb3Fnaj17BuGX2WK8x7nC1w3yTiNTPxDg2AJVEJFREcgO9AJvm1mQ4Pz/hoUZ3sfhpx2O9by3cTef3V7H/jA33bkx6uJo8GgJrRGS/iGwVkW0istWdE4rIFBwTS1URkTgRGaCq14BhwEJgJzBNVXe4c3xjXBHsfKx3fN9wzly4yj/XXmLUnB3EX7bHeo1xhasd5u0z6oSq2juF8vnYuFkmi7WrXorGFYoxfNJSPl9zkIU7fuOVLtVpZ6P1GpMql+48VPVXVf0VuIhjEqgbizFer0BgAH3D8vDd0CYUDAxg0JexDPkylt/OXvJ0aMZkW66+Yd5FRPYCB4DlwEHg+0yMy5gsV69cEeY90Yy/RVVh2e4TtH13OV+uOcj16/Y9yZikXO3zeBVoBOxR1VCgNbA206IyxkMC/P34a2RFFg5vQa2yhfi/2Tvo8fFqdv92ztOhGZOtuJo8rqrqScBPRPxUdRmQ7FuHxviCkOL5+WpAQ959oDYHT17gnrEr+feCXVy6ak9lGQOuJ48zIhKE423ur0VkDM5xrozxVSLCffXKsHhES+6tU5oPY/YTNXoFP+77w9OhGeNxriaPe3F0lj8FLAD2A50zKyhjspOi+XPzzgO1+fpRx4g5D366jhFTf+Jk/GUPR2aM57j6tNV5VU1Q1Wuq+rmqjnU2YxmTYzStWJwFw1swrFVF5mw5Spt3lzM9Ns6mvzU5kqtPW50TkT+dyyURSRCRPzM7OGOym8AAf55pX4X5TzanfIkgnvl2Cw9+uo5ffo/3dGjGZClX7zwKqGpBVS0I5AW6Ax9mamTGZGOVgwvw7eDG/KtbDbYdOUvUmJW8v2QvV67ZOFkmZ0j3HObqMIsMfOvcGG/k5yc82PAuloxoSduwYN5ZtIeOY1ey/sApT4dmTKZzaXgSEbkv0aofjsd07fVbY3BMfzuuTz161DvBi7O288Ana+hVvywjO1SjUL4AT4dnTKZwdWyrxE9WXcPxhvm9GR6NMV6sVdWSLBrRgqPaB/wAABqHSURBVNGL9zJx1QEW7zzO/3UKo0vtO2+ZeMoYX+BS8lDV/pkdiDG+IF/uXPy9YzXurXMnf5+xjSejf+K7TUf45701KFcsn6fDMybDuNpsNTa1z1X1iYwJxxjfUP3OQsz4a1O+Wvsrby3cTbvRy3mydWUebR5KgH+6uxqNyXZc/VccCNQD9jqXOkBuINa5GGOS8PcT+jUJYdGIFrSsXII3F+yi09hVxP562tOhGXPbXE0etYBIVX1fVd/HMTBiHecLg59nXnjGeL87CuXlk74RjO8bzp+XrtLj49W8MHMbZy9e9XRoxrjN1eRRBCiYaD3IWeYRIlJORGaJyCQRed5TcRiTHu2ql2LRiJb0bxLKlPWHaPPucuZtPWpvqBuv5GryeAPYLCKTReRzYBPwmjsndP7CPyEi25OUR4nIbhHZ50JCqAlMV9W/AHXdicMYTwjKk4uXOocxZ1gzShUMZNg3m+k/eQOHT13wdGjGpIurb5h/hmMe85nAd0Dj22iumgxEJS4QEX9gHNABCAN6i0iYiNQUkXlJlpI45hIZICJLcQzUaIxXqVG6ELMea8pLncLYcOAUbd9bzsfL93M1wd5QN94h1eQhIneJSCEAVf0N+BNHf0cfEcntzglVdQWQ9BXcBsA+Vf1FVa8A0cC9qrpNVTslWU4A/YGXVfVu4B534jDG0/z9hL80C2Xx0y1pWbkEb3y/i87vW4e68Q6SWnuriKwDuqnqURGpAywGXsfRgX5VVR9166QiIcA8Va3hXO8BRN04noj0BRqq6rAU9q8BjAL+AOJV9ZlkthkEDAIIDg4Oj46OdidUAOLj4wkKCnJ7f2+T0+oL2aPOm45f46udVzh9SYksm4selXOTPyDzXi7MDnXOalbn9GnVqlWsqiY78V9a73nkVdWjzp8fAiap6jsi4gf85FY0GUBVtwM90thmPDAeICIiQiMjI90+X0xMDLezv7fJafWF7FHnSGDQ5Wu8t2gPn/14gG2n/Xmpcxida92RKW+oZ4c6ZzWrc8ZJq88j8b/Yu4ElAKqa0Q2zR4CyidbLOMuMyVGC8uTi/zo5OtTvKBTIE1M20++zDRw6aR3qJntJK3ksFZFpzmlniwBLAUTkDuBKBsaxAagkIqHOvpRewJwMPL4xXuVGh/qozmFs+vU0bd9bzocx+6xD3WQbaSWP4cAMHAMhNlPVG281lQJecOeEIjIFWANUEZE4ERmgqteAYcBCYCcwTVV3uHN8Y3yFv5/wSNNQFo9oyd1VS/LvBbu5Z+xKNh60Id+N56Xa56GO3vRbeppFpJOqznP3hKraO4Xy+cB8d49rjK8qVSiQjx4KZ8nO47w0ewc9Pl5D7wZleS6qKoXzufXQozG3zZ0R2v6R4VEYY9LUulowi0a0YFCL8kzbGEfrd5Yza/MRe0PdeIQ7ycMmJjDGQ24M+T53WDPKFM3H8Kk/0Xfieg78cd7ToZkcxp3kMTjDozDGpEvYnQWZMbQJr95bnS2Hz9B+9ArGLtnL5WsJng7N5BAuJw8RaSIifYCqIvKwiDyciXEZY9Lg7yf0bRzCkqcdc6i/u2gPHcasZM3+k54OzeQALiUPEfkSeBtoBtR3Lsm+dWiMyVo35lCf3L8+VxOu03vCWp6etoVT5zPyaXpjbuXqHOYRQJhaz5wx2VZklZL8MLwl7y/dy/gVv7Bk13H+3qEa90eUsTnUTYZztdlqO453O4wx2Vje3P78Laoq859sTqWSQfztu630/GQte4+f83Roxse4mjyKAz+LyEIRmXNjyczAjDHuqxxcgKmDGvNm95rsOXGOjmNX8tbCXVy6ah3qJmO42mw1KjODMMZkPD8/oWf9crSpFsxr83cxbtl+5m45xqtda9CycglPh2e8nEvJQ1WXZ3YgxpjMUSwoD+88UJvu4aV5cdZ2+k1aT6dad9CmqI2TZdzn6tNWjURkg4jEi8gVEUkQkT8zOzhjTMZpUqE43z/ZnBFtK/PDz8cZueoiX6w5SMJ1ew7GpJ+rfR4fAL2BvUBe4FEc08YaY7xInlz+PNG6EguHtyC0kB8vzd7BfR+tZvuRs54OzXgZl18SVNV9gL+qJjjnNI9Kax9jTPYUWjw/z0YEMqZXHY6cvkCXD1bx6ryfib98zdOhGS/haof5Bec8Gz+JyL+BY7g3tIkxJpsQEe6tU5rIyiV5c+EuJq46wPxtx3i5c3XaVw+2d0NMqlxNAH2d2w4DzuOY9a97ZgVljMk6hfIF8Fq3mnw3tAmF8gYw5KtYBn6xkbjTNnuhSZlLyUNVf8Uxmu4dqvqKqo5wNmMZY3xE+F1FmPd4M17oWI0f952k7bsr+GT5fpu90CTL1aetOgM/AQuc63XsJUFjfE8ufz8GtijP4qdb0rRicV7/fhed319F7K82e6G5lavNVqOABsAZAFX9CQjNpJhuISLlRWSiiExPVJZfRD4XkQki8mBWxGFMTlK6cF4+7RfB+L7h/HnxKt0/WsPIGds4c8EGWzQOriaPq6qa9Fm+NB8OF5FJInJCRLYnKY8Skd0isk9Enk/tGKr6i6oOSFJ8HzBdVQcCXVypgDEm/dpVL8WiES15tFko0zYepvU7y5m5Oc5mLzQuJ48dzrk8/EWkkoi8D6x2Yb/JJHmkV0T8cbwj0gEIA3qLSJiI1BSReUmWkikctwxw2PmzDdZjTCbKnycXL3YKY86wppQtmo+npm7hwU/Xsf/3eE+HZjxIXPkGISL5gBeAdjg6zhcCr6rqJRf2DQHmqWoN53pjYJSqtneujwRQ1dfTOM50Ve3h/LkvcFpV54lItKr2Smb7QcAggODg4PDo6Og065mS+Ph4goKC3N7f2+S0+oLV2VXXVYk5fI1v91zhagLcUz6Ae8oHkNvfOx7rteucPq1atYpV1eTnblLVTF2AEGB7ovUewKeJ1vsCH6SyfzHgY2A/MNJZlh/4DPgIeDCtGMLDw/V2LFu27Lb29zY5rb6qVuf0Ov7nRX1iyia967l52vLfS3XFnhMZF1gmsuucPsBGTeH3aqovCab1RJWqZnp/g6qeBIYkKTsP9M/scxtjkleyQCBjetXl/vCyvDhrG30nrqdL7Tt5sVM1ShYI9HR4Jguk9YZ5Yxx9C1OAdTiarG7XERwvGd5QxllmjPEyzSoVZ8HwFnwUs5+PYvazbPcJ/hZVlT4NyuHv5x1NWcY9aXWYlwL+DtQAxgBtgT9Udbm6P0z7BqCSiIQ6hzzpBdg7I8Z4qcAAf55qW5kFw5tTq0wh/m/WdhtsMQdINXmoYxDEBaraD2gE7ANiRGSYKwcXkSnAGqCKiMSJyABVvYZjmJOFwE5gmqruuK1aGGM8rnyJIL4a0JDRPf872OI/5tpgi74qzYERRSQPcA+OIdlDgLHATFcOrqq9UyifD8x3OUpjjFcQEbrWLU2rKo7BFj9b7RhscVSXMNpXL2WDLfqQVO88ROQLHHcO9YBXVLW+qr6qqtZHYYxJUeLBFovkz82Qrzbx6OcbOXzKBlv0FWn1eTwEVAKeBFaLyJ/O5ZzNJGiMSUu9ckWYO6wpL95TjTW/nKTte8v5KMYGW/QFafV5+KlqAedSMNFSQFULZlWQxhjvlcvfj0ebl2fxiJa0qFSCNxfs4p6xK1l/wAZb9GY2oZMxJkvcWTgv4x+O4NOHIzh/OYEHPlnD36Zv4dR5G2zRG1nyMMZkqTZhwSwa0YLBLcszY9MRWr8Tw7SNh22wRS9jycMYk+Xy5c7FyA7V+M8TzalQIoi/Td9Kz0/Wsuf4OU+HZlxkycMY4zFVShVg2uDG/Lt7LfacOEfHMSt5c8EuLl6xwbKzO0sexhiP8vMTHqhflqVPR9K1bmk+itlP2/eWs3TXcU+HZlJhycMYky0UzZ+bt++vzdRBjQgM8Ocvkzcy5MtYjp296OnQTDIseRhjspWG5Ysx/4nmPNu+Cst2n6DNO8v5dOUvXLN3Q7IVSx7GmGwndy4/HmtVkUVPtaR+aFH++Z+ddP7gRzYfOu3p0IyTJQ9jTLZVrlg+PnukPh89WI/T569w30ereWHmNs5euOrp0HI8Sx7GmGxNROhQ8w4WP92S/k1CmbL+EK3fjWHW5iP2bogHWfIwxniFoDy5eKlzGHOGNaN0kXwMn/oTfSasY//v8Z4OLUey5GGM8So1ShdixtAmvNq1BtuPnqXD6JW8u2gPl67auyFZyZKHMcbr+PsJfRvdxdKnI+lYsxRjl+yl/egVrNjzu6dDyzG8InmISHkRmSgi0xOVdRWRCSIyVUTaeTI+Y4xnlCiQh9G96vL1ow3xF+HhSesZ9s0mjv95ydOh+bxMTx4iMklETojI9iTlUSKyW0T2icjzqR1DVX9R1QFJymap6kBgCNAz4yM3xniLphWL8/3w5oxoW5kffj5Om3eW8/nqgyRctw71zJIVdx6TgajEBSLiD4wDOgBhQG8RCRORmiIyL8lSMo3jv+g8ljEmB8uTy58nWlfih+EtqFOuMC/P2UHXcT+yNe6Mp0PzSZIVj7qJSAgwT1VrONcbA6NUtb1zfSSAqr6exnGmq2oP588CvAEsUtXFyWw7CBgEEBwcHB4dHe12/PHx8QQFBbm9v7fJafUFq7OvUVXW/5bAN7uu8Odl5e5yueheKTfXL5/32Tqn5Hauc6tWrWJVNSK5z3LdVlTuKw0cTrQeBzRMaWMRKQb8C6grIiOdSeZxoA1QSEQqqurHifdR1fHAeICIiAiNjIx0O9iYmBhuZ39vk9PqC1ZnX9QKGHrpKu/+sIcv1hxk62l/upcP5Ll7WuL47pkzZNZ19lTySBdVPYmjbyNx2VhgrGciMsZ4g4KBAYzqUp3u9crwwqxtfLzlLNvPr+fVrjUILZ7f0+F5NU89bXUEKJtovYyzzBhjMlzNMoWY+demPFQtN1sOn6H96BWMXmzvhtwOTyWPDUAlEQkVkdxAL2COh2IxxuQA/n5Cm7sCWPJ0S9pXL8XoxXvpMGYlq/b+4enQvFJWPKo7BVgDVBGROBEZoKrXgGHAQmAnME1Vd2R2LMYYU7JgIO/3rsuXAxqgqjw0cR1PTNnMiXP2bkh6ZHqfh6r2TqF8PjA/s89vjDHJaV6pBAuGt+DDmP18HLOfZbtO8GxUFR5seBf+fjmnQ91dXvGGuTHGZIbAAH9GtK3MguHNqV22MC/N3kG3D39kW9xZT4eW7VnyMMbkeOVLBPHlgAaM6VWHY2cvce+4Vbw8ezt/XrJ5Q1JiycMYY3DMG3JvndIsebolfRvdxRdrf6X1O8uZs+WozRuSDEsexhiTSMHAAF65twazH2vKHYUCeWLKZh6etJ6Df5z3dGjZiiUPY4xJRq0yhZn516a80qU6mw+dod3oFYxZvJfL1+zdELDkYYwxKfL3E/o1CWHJ0y1pFxbMe4v3EDXa3g0BSx7GGJOm4IKBfNCnHl/8pQHXE78bkoPnDbHkYYwxLmpRuQQLh7fgydaVWLD9N1q/s5zJPx7IkfOGWPIwxph0CAzw56m2lVn4lGPekFFzf+becavYcjhnzRtiycMYY9wQWjw/X/ylAR/0qcuJPy/T9cMfeXHWNs5ezBnvhljyMMYYN4kInWrdyZKnW/JIkxC+WXeI1u/EMHNznM+/G2LJwxhjblOBwABe7lydOcOaUaZIPp6auoXeE9ay9/g5T4eWaSx5GGNMBqlRuhAzhjbhtW412XnsHB3GrOT1+Ts5f/map0PLcJY8jDEmA/n5CX0almPp0y25r15pPlnxC63fWc732475VFOWJQ9jjMkExYLy8O8etfluaBOK5M/N0K830X/yBg6dvODp0DKEJQ9jjMlE4XcVYe6wpvxfpzA2HDhF2/eW8+4Pu7lwxbubsrJ98hCR8iIyUUSmJynPLyIbRaSTp2IzxhhX5PL3Y0CzUJY8HUn76qUYu3Qfd7+9nJmb47jupS8YZmryEJFJInJCRLYnKY8Skd0isk9Enk/tGKr6i6oOSOaj54BpGRmvMcZkplKFAhnbuy7fDW1MyYJ5eGrqFu77aDWxv572dGjpltl3HpOBqMQFIuIPjAM6AGFAbxEJE5GaIjIvyVIyuYOKSFvgZ+BE5oZvjDEZL/yuosz6a1Pevr82R89cpPtHq3l8ymbiTntPf4hkdu+/iIQA81S1hnO9MTBKVds710cCqOrraRxnuqr2cP78LyA/juRzEeimqteTbD8IGAQQHBwcHh0d7XYd4uPjCQoKcnt/b5PT6gtW55wiO9b50jVl/oGrfH/A8WZ6+5AA7ikfQN5cGTOP+u3UuVWrVrGqGpHcZ7luKyr3lAYOJ1qPAxqmtLGIFAP+BdQVkZGq+rqqvuD87BHgj6SJA0BVxwPjASIiIjQyMtLtgGNiYrid/b1NTqsvWJ1ziuxa5yjg6JmLvLVwNzM3H2HtCWFIywr0alCOoDy392s6s+rsieSRLqp6EhiSwmeTszYaY4zJHHcWzst7PevwSJMQ3vh+F//8z07eX7qPfo3v4qFGd1GyYKCnQ7yFJ5LHEaBsovUyzjJjjMnxapctzJRBjdh86DQfL9/P2KX7GBezn1ZVStKzflkiq5QgwN/zD8p6InlsACqJSCiOpNEL6OOBOIwxJtuqW64In/SN4MAf55m64TDfbYpj8c7jFM4XQFT1UnSqdScNQouSO5dnEkmmJg8RmQJEAsVFJA54WVUnisgwYCHgD0xS1R2ZGYcxxnir0OL5eb5DVZ5uV5nlu39n7tajzN1ylOgNh8nt70fVOwpQo3QhKpUMolzRfJQtmo9ShQIpkCcXIhnT6Z6cTE0eqto7hfL5wPzMPLcxxviSAH8/2oQF0yYsmEtXE1ix53diD51mW9xZ5m45yrlLt76xnjuXH8Xz56ZqoWtkxjMC2b7D3BhjzK0CA/xpV70U7aqXAkBVOXn+CodOXeDwqQuc+PMyf8Rf5vf4y/jHZ87rcJY8jDHGy4kIxYPyUDwoD/XKFbnls5iYmEw5p+e77I0xxngdSx7GGGPSzZKHMcaYdLPkYYwxJt0seRhjjEk3Sx7GGGPSzZKHMcaYdLPkYYwxJt0yfTKo7EBEfgd+vY1DFAf+yKBwvEFOqy9YnXMKq3P63KWqJZL7IEckj9slIhtTmk3LF+W0+oLVOaewOmcca7YyxhiTbpY8jDHGpJslD9eM93QAWSyn1ReszjmF1TmDWJ+HMcaYdLM7D2OMMelmycMYY0y6WfJIhYhEichuEdknIs97Op7MICJlRWSZiPwsIjtE5ElneVERWSQie51/FknrWN5GRPxFZLOIzHOuh4rIOuf1nioiuT0dY0YSkcIiMl1EdonIThFp7MvXWUSecv6b3i4iU0Qk0BevsYhMEpETIrI9UVmy11Ucxjrrv1VE6rl7XkseKRARf2Ac0AEIA3qLSJhno8oU14CnVTUMaAQ85qzn88ASVa0ELHGu+5ongZ2J1t8E3lPVisBpYIBHoso8Y4AFqloVqI2j7j55nUWkNPAEEKGqNQB/oBe+eY0nA1FJylK6rh2ASs5lEPCRuye15JGyBsA+Vf1FVa8A0cC9Ho4pw6nqMVXd5Pz5HI5fKKVx1PVz52afA109E2HmEJEywD3Ap851Ae4Gpjs38ak6i0ghoAUwEUBVr6jqGXz7OucC8opILiAfcAwfvMaqugI4laQ4pet6L/CFOqwFCovIHe6c15JHykoDhxOtxznLfJaIhAB1gXVAsKoec370GxDsobAyy2jgb8B153ox4IyqXnOu+9r1DgV+Bz5zNtV9KiL58dHrrKpHgLeBQziSxlkgFt++xomldF0z7PeaJQ8DgIgEAd8Bw1X1z8SfqeN5bp95pltEOgEnVDXW07FkoVxAPeAjVa0LnCdJE5UvXWdnG/+9OJLmnUB+/rdpJ0fIrOtqySNlR4CyidbLOMt8jogE4EgcX6vqDGfx8Ru3s84/T3gqvkzQFOgiIgdxNEfejaM/oLCziQN873rHAXGqus65Ph1HMvHV69wGOKCqv6vqVWAGjuvuy9c4sZSua4b9XrPkkbINQCXn0xm5cXS2zfFwTBnO2dY/Edipqu8m+mgO0M/5cz9gdlbHlllUdaSqllHVEBzXdamqPggsA3o4N/O1Ov8GHBaRKs6i1sDP+O51PgQ0EpF8zn/jN+rrs9c4iZSu6xzgYedTV42As4mat9LF3jBPhYh0xNE27g9MUtV/eTikDCcizYCVwDb+2/7/dxz9HtOAcjiGs39AVZN2ynk9EYkEnlHVTiJSHsedSFFgM/CQql72ZHwZSUTq4HhAIDfwC9AfxxdIn7zOIvIK0BPHE4WbgUdxtO/71DUWkSlAJI6h148DLwOzSOa6OhPpBzia8C4A/VV1o1vnteRhjDEmvazZyhhjTLpZ8jDGGJNuljyMMcakmyUPY4wx6WbJwxhjTLpZ8jAmnURktfPPEBHpk8HH/nty5zImu7FHdY1xU+J3RNKxT65EYysl93m8qgZlRHzGZCa78zAmnUQk3vnjG0BzEfnJOXeEv4i8JSIbnHMlDHZuHykiK0VkDo63nBGRWSIS65xvYpCz7A0co8D+JCJfJz6X843gt5xzU2wTkZ6Jjh2TaJ6Or50vghmTqXKlvYkxJgXPk+jOw5kEzqpqfRHJA/woIj84t60H1FDVA871vzjf+M0LbBCR71T1eREZpqp1kjnXfUAdHPNwFHfus8L5WV2gOnAU+BHHGE6rMr66xvyX3XkYk3Ha4Rg36Cccw7sUwzHpDsD6RIkD4AkR2QKsxTFQXSVS1wyYoqoJqnocWA7UT3TsOFW9DvwEhGRIbYxJhd15GJNxBHhcVRfeUujoGzmfZL0N0FhVL4hIDBB4G+dNPDZTAvb/2mQBu/Mwxn3ngAKJ1hcCQ51D3CMilZ0TLiVVCDjtTBxVcUz/e8PVG/snsRLo6exXKYFjVsD1GVILY9xg31CMcd9WIMHZ/DQZx5wgIcAmZ6f17yQ/zekCYIiI7AR242i6umE8sFVENjmHib9hJtAY2IJjYp+/qepvzuRjTJazR3WNMcakmzVbGWOMSTdLHsYYY9LNkocxxph0s+RhjDEm3Sx5GGOMSTdLHsYYY9LNkocxxph0+39i+gxJTehDEAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import scipy.io as sio\n",
    "import os\n",
    "%matplotlib inline\n",
    "mse_records_dict = sio.loadmat('results/CentralizedADMM.mat')\n",
    "\n",
    "plt.semilogy(mse_records_dict[\"mse\"][0], label=\"Centralized ADMM\")\n",
    "plt.title(\"Centralized ADMM\")\n",
    "plt.xlabel(\"iteration\")\n",
    "plt.ylabel(\"Mean-Square-Error\")\n",
    "plt.grid(\"on\")\n",
    "plt.legend()\n",
    "dirname = 'images'\n",
    "if not os.path.exists(dirname):\n",
    "    os.makedirs(dirname)\n",
    "plt.savefig(os.path.join(dirname, 'centralized_admm.png'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "You can find the full code in the `CentralizedADMM.py`\n",
    "\n",
    "**Exercise**: Modify the loss function, (sub-)gradient, and proximal step for $\\ell_1 + \\ell_2$ cases."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}